function IRF = zFC_IRFs(Ahat, P, T_irf, shocked_variab, shock_size)


% -----------------------------------------------------%
% DESCRIPTION OF THE FUNCTION
% this function computes impulse response functions in a VAR 

% Generates a k x T_irf matrix containing the impulse response functions of
% the k variables to a shock of value "shock_size" to variable "shocked_variable". 
% Ahat contains the coefficients of the autoregressive VAR model, P the
% matrix that maps structural shocks into reduced form shocks
% -----------------------------------------------------%

    assert( size(Ahat,1) <= size(Ahat,2) )

    k = size(Ahat,1); % number of variables in the underlying VAR model
    p = size(Ahat,2)/k; % number of lags in the underlying VAR model

    Spseudo = zeros(k,T_irf);
    Spseudo(shocked_variab,1) = shock_size;
    Rpseudo = P*Spseudo;

    IRFestim = zeros(k,T_irf); 
    IRFestim(:,1) = Rpseudo(:,1); 
    for j=1+1:T_irf
        step1 = [fliplr(IRFestim), zeros(k,p-1)];
        step2 = reshape(step1,k*T_irf+k*(p-1),1);
        step3 = step2(k*(T_irf-j+1)+1:k*(T_irf-j+1+p),1);
        IRFestim(:,j) =  Ahat*step3 + Rpseudo(:,j);  
    end % do it better

    IRF = IRFestim;

end